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I. INTRODUCTION 

Nuclear theory seeks to formulate a consistent frame- 
work for nuclear structure and reactions. Nuclear density 
functional theory (DFT) is a particularly useful tool for 
describing the ground-state properties of nuclei across the 
nuclear chart 1-4 . This approach is based on the pio- 
neering works of Skyrme [5 and Vautherin [SI [7] and 
might theoretically be based on the theorems by Ho- 
henberg and Kohn |5]. It has been applied through 
the self-consistent mean-field computations with density- 
dependent energy functionals Formally, the energy- 
density functional results from the Legendre transform 
of the ground-state energy as a functional of external 
one-body potentials [TU] . While DFT allows for a simple 
solution to the quantum many-body problem, the con- 
struction of the functional itself poses a challenge [TTlll2j . 
Various efforts aim at constraining the functional from 
microscopic interactions |13j or devising a systematic ap- 
proach |14) . In this paper we explore the development 
of an energy functional that employs shell-model occu- 
pations instead of the nuclear density. The minimization 
of this functional is technically simpler than DFT, and 
global mass-table fits of the functional require only mod- 
est computational resources. 

In nuclear physics, energy-density functional theory is 
a practical tool that is popular because of its compu- 
tational simplicity and success [2 IE] ■ The universality 
of the functional (i.e., the possibility to study nucleons 
in external potentials) is seldom used [TB]. This dis- 
tinguishes nuclear DFT from DFT in Coulomb systems 
and makes alternative formulations worth studying. For 
computational simplicity, we would like to maintain the 
framework of an energy functional. However, there is no 
need to focus only on functionals of the density. The 
description and interpretation of nuclear systems are of- 
ten based on shell-model occupation numbers rather than 
densities [T7]. For instance, in nuclear structure one is 
often as much interested in the occupation of a given 
shell-model orbital or the isospin-dependence of the effec- 
tive single-particle energies as in the shape of the density 
distribution. This approach is natural for shell-model 



Hamiltonians that are based on single-particle orbitals. 
The purpose of this paper is to develop a nuclear energy 
functional constructed in a shell-model framework. We 
show that an occupation number-based functional can 
be employed in global mass-table computations and can 
perform at a competitive level, with a reasonable number 
of parameters and relatively short computation time. 

Many approaches to the nuclear energy-density func- 
tional are empirical [ISHST] . but guidance can be found 
from analytical solutions to simpler models. The com- 
plete ground-state energy of a quantum many-body sys- 
tem as a functional of any external potential is available 
for only a few solvable or weakly interacting systems. 
Furnstahl and colleagues P^H^ derived energy-density 
functionals for dilute Fermi gases with short-ranged in- 
teractions. For Fermi gases in the unitary regime, simple 
scaling arguments suggest the form of the energy-density 
functional [5SH55] . We similarly use scaling arguments in 
the development of our functional. 

In recent years, significant progress has been made 
in developing global nuclear energy-density functionals; 
see, for example, Refs. [251 133 • The Skyrme-Hartree- 
Fock-Bogoliubov mass functionals by Goriely et al. [2^ 
achieve a least-squares error of about x — 0.58 MeV 
with 14 parameters fit to nuclei with N,Z > 8. The 
12-parameter UNEDF functionals UNEDFnb and UN- 
EDFO of Ref. [30_ have a least-squares deviation of 
X — 0.97 MeV and x = 1-46 MeV, respectively, when 
fit to a set of 72 even-even nuclei. The fitted function- 
als, when applied to 520 even-even nuclei of a mass ta- 
ble, yield least-squares deviations of x = 1-45 MeV and 
X = 1.61 MeV, respectively. The refit [311 of the Skyrme 
functional Sly4 to even-even nuclei exhibits a root-mean- 
square deviation of about x = 1-7 MeV. For a recent 
review on this matter, we refer the reader to Ref. |15j . 
The work presented here provides an approach based 
on Hohenberg-Kohn DFT, utilizing degrees of freedom 
that are natural to the nuclear shell-model, and gives a 
global description of the nuclear chart for all nuclei with 
N,Z>8. 

We emphasize the difference between an energy 
functional for nuclear masses and a mass formula. 
Mass formulae, such as the finite-range droplet model 
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(FRDM) and Duflo-Zuker mass formula [33], calcu- 
late nuclear masses from an assumed nuclear structure. 
The mass formula by Duflo-Zuker, for instance, computes 
the binding energy in terms of the shell-model occupa- 
tions. The occupations are parameters — not variables — 
and their values are taken from the noninteracting shell- 
model (with modifications due to deformations). The 
form of the Duflo-Zuker mass formula is motivated by 
semiclassical scaling arguments, which ensure that the 
binding energy scales linearly with the mass number in 
leading order; see Ref. [33] for details. The Duflo-Zuker 
mass formula achieves an impressive root-mean-square 
deviation of x = 0.35 MeV with 28 fit parameters. In an 
occupation number-based energy functional, the ground- 
state energy results from a minimization of a functional, 
and the occupation numbers are variables that minimize 
the functional. In other words, the structure of a nu- 
cleus results from the minimization of a functional and 
is not assumed. In our phenomenological construction of 
an occupation number-based energy functional, we em- 
ploy several scaling arguments and terms from the Duflo- 
Zuker mass formula. Recently, the energy functional (in 
terms of occupation numbers) was constructed for the 
pairing Hamiltonian [35H37j and the three-level Lipkin 
model [381 ES] ■ The analytical insights gained from these 
simple models also enter the construction of terms for 
our occupation number-based functional. 

This paper is organized as follows. In Sect. |ll] we in- 
troduce the form of the energy functional and discuss the 
relevant terms. In Sects. [TTC] and [Till we provide a statis- 
tical analysis of the functional and a description of the 
minimization routines used. We study the performance 
of the functional in calculating nuclear ground-state en- 



ergies and radii in Sect. IV We conclude with a summary 
in Sect. |Vl 



II. OCCUPATION NUMBER-BASED ENERGY 
FUNCTIONAL 

In this section, we introduce the theoretical foundation 
for functionals based on occupation numbers and discuss 
the terms that enter the functional. 



Current nuclear DFT models, such as HFB-17 [25], are 
based on Skyrme forces and employ densities p(r) and 
currents depending on spatial coordinates. The use of 
such objects seems a natural choice, particularly in quan- 
tum chemistry and condensed matter systems where one 
is interested in the electronic structure in the presence 
of ions. For atomic nuclei, however, other representa- 
tions might also be interesting or more natural, and we 
consider functionals based on occupations of shell-model 
orbitals. 

To present our concept formally, we consider the 
ground-state energy E{ea) of an A-body system de- 
scribed by a shell-model Hamiltonian with single-particle 
energies Sa, a — 1, 2, 3, . . ., 



(2) 



Here, a|j creates a fermion in the shell-model orbital a, 
and V is the nuclear interaction. The ground-state en- 
ergy is a (complicated) function of the single-particle en- 
ergies, and the occupation number of the orbital with 
label /3 is 



_d_ 



E{e^). 



(3) 



A Legendre transformation of the function E{ea) yields 
the occupation number-based energy functional [lOj 



(4) 



From a technical point of view, the one-body potential — 
or the mean-field — has been employed as an external po- 
tential in DFT. Thus, there is a formal analogy between 
Hohenberg-Kohn DFT and occupation number function- 
als. The Hohenberg-Kohn theorems remain valid when 
we exchange w(r) o £a and p(r) o ria- However, the 
explicit construction of the occupation number-based en- 
ergy functional as a Legendre transform is possible only 
for exactly solvable systems [35] [3S] . In the following sec- 
tion, we employ primarily phenomenological arguments 
for the terms that enter the functional. 



Theoretical Basis 



B. Form of the Functional 



We first clarify our work within Hohenberg-Kohn 
DFT ^8. • According to Hohenberg and Kohn, for a many- 
body system with density p in an external potential v 
there exists a functional of the density F-uKip) indepen- 
dent of V such that the ground-state energy is 



(1) 



Here the minimization is performed over all one-body 
densities p{r). 



The occupation number-based functional is guided by 
semiclassical scaling arguments that guarantee nuclear 
saturation |331l34j . In the harmonic oscillator basis, a nu- 
cleus of A nucleons has about A^^^ occupied shells below 
the Fermi level, and about A"^/"^ nucleons are in the high- 
est energetic, completely filled shell. We are interested 
in expressing the functional as a sum of simple functions 
that depend on occupation numbers. Consider, for exam- 
ple, the energy from a harmonic mean-field potential for 
the protons. Let Zp denote the proton occupation num- 
ber of the pth oscillator shell. The maximum occupation 
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of the pth shell is 



pzp^p , 



(5) 



where ^ indicates "scales as." We use semiclassical ar- 
guments to enforce saturation. The total energy of the 
harmonic mean-field scales as 



E 



^1/3 

pzp ~ I dp 





(6) 
(7) 



Thus, the multiplication of such an energy term with a 
factor ofA~^/^ ensures saturation, that is, a scaling of the 
energy with A in leading order. The simple arguments 
based on filling the major shells of the spherical harmonic 
oscillator yield incorrect shell closures. As a remedy, we 
follow Duflo and Zuker 33 and consider shells where 
the high-j intruder orbital of the shell p + 1 enters the 
shell p. Since a single j orbital of the shell p holds on 
the order of p nucleons — while a major shell holds about 
p^ nucleons — our scaling arguments are unchanged. In 
the modified shell-model, the resulting shell closures have 
occupation numbers 



2,6,14,28,50,82,126,184,. 



(8) 



In what follows, we employ a model space consisting 
of 15 major shells. The shell p has a dimension dp — 
[p -I- l)(p + 2), (i.e., it can hold up to dp protons and 
dp neutrons). We employ suitable factors of dp and p to 
ensure saturation. 

We write the functional as a sum of macroscopic 
Coulomb and pairing terms plus a microscopic interac- 
tion z) that depends on the occupations Zp of proton 



shells and Up of neutron shells 



F{c; n, z) = c, 



Z{Z-l) 

^1/3 



Cp 



+ J"(c;n,z). (9) 



Here, 5 is 1, 0, and -1 for even-even, odd mass, and odd- 
odd nuclei, respectively, and Z is the charge number. The 
17 fit parameters of the full functional are denoted as the 
shorthand coefficients 



C = {Ci, C2, . . . , Cii,Cc, Cp, Cs, Cas, Css, C^}, (10) 



while 



Z = {Zl, Z2, . . .} 

n = {ni,n2, . . .} 



(11) 
(12) 



denote the occupations of proton and neutron shells, re- 
spectively. The microscopic functional J-{c] n, z) has the 
form 



J-{c] n, z) — huiV + T]. 



'kin 



l2B 



D 



+DiB +T+ Af4Bcx + ival 

+hCb{D2B+ L + L). 



(13) 



Below, we explain the individual terms entering the func- 
tional. The functions Hlu and huj are 



huj{Cs,Cas,Css) = i-CsA 



T(T+ 1) 



1 + CssA-^/3 A^ 

huj{Cs) = l+C^A^^/^, 



, (14) 
(15) 



and they account for a smooth dependence on the mass 
and the isospin T — N — Z (with N being the neutron 
number). Each of the terms of the functional (131 scales 
at most as A. Thus, the functions (14 1 and (151 include 
smooth surface corrections. The last term of Eq. (14) 
provides both a volume- and surface-symmetry correc- 
tion [40, 41J. 

We now discuss the individual terms of the func- 



tional (13). We begin with the volume terms that scale 



as A. These are 



V 

Tkin 
-^2B 



ciA, 

C2A-^/^^p{zp + np), 



(16) 
(17) 



C3 



Zp(Zp 1) 



p 



np{np - 1) 2 

p p 



(18) 



D 



p 

E 



j3/2 



{zp - dp/2f 



d] 



jj^{n,^d,/2)' 



4B 



ceE 



q \ -9 



dl 



(19) 
(20) 



^ = ^^E v(^p + £)(^*p + i) 



+ '^{np + e){zp + l) 



T 



cg E V " ^p)^ + ' 
p 



(21) 
(22) 



Note that e — 10 ^ is a regularization parameter in 



Eqs. (16)-(22), to ensure differentiability. The term V of 



Eq. ( 16 ) is a smooth volume term and accounts for the 



bulk energy. The term Tkin of Eq. ( 17 ) contains a contri- 



bution of the single-particle kinetic energies. This term 
exhibits shell effects and thereby corrects the smooth 
behavior of the volume term V . Effects of two-body 
interactions are captured by the terms beginning with 
Eq. (18). The two-body term /2B is motivated by the 



monopolc Hamiltonian in the mass formula by Duflo and 
Zuker [33]. Here, the fermionic nature in the proton- 
proton and neutron-neutron interaction is evident. 
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The term D of Eq. (19) is technically a two-, three- 



and four-body interaction and serves to include deforma- 
tion effects. The form of this term has a linear onset 
for almost-empty and almost-filled shells and assumes 
its maximum at mid-shell with half-filled occupations. 
In practice, this term behaves similar to the function 
min(np,dp — Hp) but does not suffer from discontinu- 
ity and lack of differentiability. Similarly, the four-body 
term D^b of Eq. ( |20| ) also accounts for deformation ef- 
fects. 



The term of Eq. (21 1 is motivated by the analytical 
results in Refs. [331 l38j. For a two- level system with 
interactions between the levels, the minimum energy re- 
sults in a fractional occupation of the higher level due to 
a square root singularity in the functional. Such square 
root terms are present in the functional of the three-level 
Lipkin model [3£j and the pairing Hamiltonian [35j |36] . 
The inclusion of these analytically motivated terms con- 
siderably improves the functional's reproduction of ex- 
perimental energies. 



The term T of Eq. ( 22 ) is an isospin term and includes 



the isospin contributions of individual orbitals. Not sur- 



prisingly, Eq. ( 22 ) turns out to be relevant for the ac- 
curate descriptions of nuclei far from the N = Z line. 



Zp\ » £ we 



For computational purposes Eq. ( 22 1 is written in square 
roots rather than absolute values. For \n. 
have ^y {up — Zp)'^ + e 
practically a differentiable approximation of the absolute 
value. 



In-n — zJ, and this term is thus 



The volume terms of Eqs. ( 16 ) to ( 22 ) also have surface 
corrections due to the discrete sums. However, we also 
need to include surface corrections independent of the 
volume terms. Technically, surface terms scale as A^^^. 
We employ 



Zp[dp Zp) ^ Up^dp Up) \ (23) 



j3/2 



j3/2 



D2B = 

p \ "-p ^P / 

L = 08^-2/3 ^ {zp^Zp + 1 + np^Up + 1) (24) 



4Bc 



ClO 



n23 



{zp{zp~l){zp~2){zp~?,)f 



+ {np{np - l){np - 2){np - 3))"^ 



Cll 



{df - n/)n/+i + I 



+ ^{df - zf)zf+i+e^y 



(25) 



(26) 



The term Z?2b of Eq. ( 23 1 accounts for microscopic pair- 
ing contributions. Its form is motivated by Talmi's se- 
niority model |42j . again with a maximum contribution 
at mid-shell as with the volume term of Eq. (19 1. Like- 



wise, the term L of Eq. (24 1 is a surface correction term 
to the volume term (21). 



The term in Eq. ( 25 ) accounts for higher-order effects 



and resulted from an efhcient scheme to identify relevant 



contributions that we describe in Sect. II D The powers 
of p ensure the desired A^^^ scaling. 

The term Lvai of Eq. ( 26 1 is unique as it involves only 



the highest occupied shell(p — /) in the noninteracting 
shell-model; that is, it treats the shell immediately below 
(with occupation Uf) and above (with occupation ?i/+i) 
the Fermi surface. In the naive shell-model, these shells 
contain about A^^^ nucleons, and Lvai thus is a surface 
term. This term accounts for energy from particle-hole 
excitations and turns out to be useful. We note that the 



forms of Eqs. (16l-(26l are symmetric about exchange 



of protons and neutrons, Zp o Up. Therefore, the only 
isospin breaking term of our functional is the Coulomb 
term of Eq. 



C. Numerical Implementation 

We consider one of the A^pts nuclei of the nuclear chart 
and label it by i. This nucleus has neutron number Ni 
and proton number Zi. The minimization of the func- 
tional ([9| with respect to the occupations n, z yields the 
ground-state energy 



Ei{c; n*{c),z*{c)) 



min F{c;n,z). (27) 



We refer to this minimization as the lower-level optimiza- 
tion. In Eq. (27), n*{c),z*{c) denote the optimal occupa- 



tion numbers (we suppress that they depend on the label 
i) and are taken from the domain 



p 



{n,z) : 



N,: 



(28) 



p 



Zi'jO < Up, Zp < dp 



From a mathematical point of view, we deal with an 
optimization problem that includes linear equality con- 
straints and bound constraints. The functional ^ is 
nonlinear but twice continuously differentiable in the oc- 
cupation numbers n, z over the domain V. Furthermore, 
given the coefficients c, we have algebraic derivatives of 
F with respect to n, z and can use these to solve each 
lower-level problem. We used the sequential quadratic 
programming (SQP) routine FFSQP [33] to solve the con- 
strained optimization problem (27). 



We provided FFSQP with analytical expressions for the 
gradient of F with respect to n and z. These gradients 
are used by FFSQP to build second-order approximations 
of the objective while remaining within the constrained 
domain Vi. The value of F is iteratively reduced until 
further local decreases would be infeasible with respect to 
Vi. This provides a local minimum in the neighborhood 
of our starting point at the naive level-filling. We found 
that the additional effort of employing the derivatives 
with respect to the occupation numbers provided signif- 
icant advantages, in terms of speed and stability, over 
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other optimization routines such as COBYLA 
not require derivatives to be available. 



that do 



As a starting point, we provide FFSQP with the occu- 
pation numbers n and z corresponding to the naive level 
filling, which are feasible for the domain Pj. In practice, 
a local solution is found within about a dozen function 
and gradient calls. 

The fit coefficients c are then determined by minimiz- 
ing the sum of squared residuals 



iVnts — 



(29) 



pts . T 



as a function of c. Here, Ef^^ denotes the measured 
ground-state energy of the nucleus i. Technically, we deal 
with a multi-level optimization problem because each of 
the theoretical energies Ei is a solution of a lower-level 
optimization of the functional ^ over the occupation 
numbers n and z. Obtaining the occupations from an op- 
timization is in contrast to the mass formula |33| . where 
occupations are fixed and input by hand. We refer to 
the optimization of the coefficients c as the upper-level 
optimization problem. This problem is solved with the 
program POUNDerS (Practical Optimization Using No 
Der ivatives for sums of Squares), and details are pre- 
sented in Sect. IIIII 



Note the computational complexity present in the min- 



imization of the chi-square ( 29 1 . For each of the A'pts 



nuclei we must perform a lower-level minimization of the 
functional to determine occupation numbers that mini- 
mize the functional and yield the ground-state energy Ei. 
The occupation numbers resulting from this lower-level 
optimization can be very sensitive to the coefficients c 
of the functional, which in turn are determined in the 
upper-level chi-square minimization. To reduce the wall 
clock time of each evaluation, we can perform the 
lower-level minimizations in parallel, using as many as 
A'pts processors. Further details of the minimization 



are presented in Sect. Ill 



The computation of ground-state energies via the func- 
tional ^ with coupled minimizations differs consider- 
ably from the uncoupled case. If the proton and neu- 
tron occupations are kept fixed to the naive level filling, 
Eq. ([9 ) becomes a mass formula. Performing a chi-square 
fit with this fixed filling yields coefficients c that can then 
be used when minimizing the functional with respect to 
occupations (that is, the lower-level minimization is done 
independently of the upper-level minimization). The re- 
sulting error, x — 3-38 MeV, is one order of magnitude 
larger than the residual error of the mass formula by Du- 
flo and Zuker [33]. When including the lower-level op- 
timization in the chi-square fit (and hence dealing with 
a multi-level optimization), we obtain a considerably re- 
duced X = 1-31 MeV. 



D. Determination of Relevant Terms 

Unfortunately, there is no recipe available for the con- 
struction of the occupation-number-based energy func- 
tional, and the reader may wonder how we arrived at the 
particular form ( 13 ) of the functional. As described in 



Sect. |II B[ the guiding principles are saturation, mean- 
field arguments, and insights from analytically solvable 
systems. However, these arguments do not fully con- 
strain the functional, and we need a more systematic ap- 
proach to identify terms that should enter into the func- 
tional. 

Bertsch et al. |3TI employed the singular value decom- 
position and studied the statistical importance of various 
linear combinations of terms that enter the functional. 
This method identifies the relative importance of pos- 
sible combinations of terms and truncates search direc- 
tions that are flat in the parameter space. Along these 
ideas, we employ a method that chooses new functional 
terms based on their correlation to terms already present. 
New terms are chosen to provide a relatively independent 
"search direction" in the parameter space of the coeffi- 
cients c in which the of Eq. (29) is minimized. This 



approach is presented in detail in Subsect. |II D 1[ 

In mass formulae, the addition of new terms (and new 
fitting parameters) yields a chi-square that is a monotone 
decreasing function with the number of fit parameters. 
For a functional, however, matters are different. Here, 
the addition of a new term to the functional guarantees 
a lowered chi-square only if the term has a perturbatively 
small effect. This point is discussed in Subsect. |IID 2| 



1. Correlation Test 

We first describe our method for selecting new terms to 
be included in the functional, beyond those strongly mo- 
tivated by scaling arguments, mean-field arguments, and 
solutions to simple analytic systems. We seek a system- 
atic method for determining new terms that will provide 
further insight into the physical system and decrease the 
overall chi-square. 

Consider the addition of a term c//, with new fit pa- 
rameter c/, to the functional 



M 



Z(Z 1^ s 



Faic;n,z 
(Z 



This functional contains M terms depending on occupa- 
tion numbers, denoting the occupation numbers 
that minimize for a given nucleus. The occupation 
numbers depend on the nucleus i under consideration, 
but we suppress this dependency. One can expect that 
the addition of the term c// to Fq will be useful in low- 
ering the chi-square only when it is independent of the 
M terms already included in the functional. 
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For identification of a new search direction we compute 
the correlation coefficient 



obtained from the functional Fq. Let us assume that 



cov(/a,/) 



Here the covariance is 

cov(/„,/) = (/„/) -(/„)(/), 



(31) 



(32) 



and the average (•) is computed with respect to the A'pts 
nuclei of the nuclear chart. Similarly, the standard devi- 
ations are 



3/ 



ViP) - if)' and Sf^ = J if J) - {f^y. (33) 



In the computation of the averages, the terms / and 
fa are evaluated at the occupations (n*(c), 2;*(c)) that 
minimize Eq. (30 1. Here (n*(c), z*(c)) depend on the 
coefficients c that minimize the based on the ener- 
gies in Eq. (30 1. Wc note that the correlation coeffi- 



cient is independent of the size of the coefficient c/ of the 
new term under consideration, though dependent on the 
other fit coefficients c through the optimal occupations 
(n*(c), 2;*(c)). Should the correlation be sufficiently low 
for all included terms, the new term c// becomes a can- 
didate and is tested for its performance in lowering the 
X value. 

This approach allows us to probe many different forms 
of functional terms and then scan through several hun- 
dred iterations without the time-consuming and compu- 
tationally expensive aspects of performing a full mini- 
mization of the chi-square for each new term under ques- 
tion. In this way, we systematically grow the functional 
term by term. We started from an initial base of about 
350 different terms and found that only 18 of them were 
weakly correlated to the existing terms and had the po- 
tential to significantly decrease the least-squares error. 
Of these 18 terms, 15 were seen to be simply higher-order 
corrections of three primary forms. We further reduced 
the set of possible terms through physical arguments and 
preliminary fits. This approach showed Eq. (22) to be 



the best choice for lowering our least-squares deviation. 
Equations ( 25 1 and ( 26 1 were determined similarly, from 
very large sets of possible terms. In this way we suc- 
cessfully lowered the functional's least-squares error to a 
meaningful x = 1-31 McV with 17 fit parameters. 



2. Perturbative Test 



Assume again that our functional Fq is as in Eq. (301, 



and that we consider the addition of a new term c// (with 
the new fit coefficient c/ taken to be the mean value of 
the currently determined fit coefficients, where there are 
no statistical outliers). We will consider the general case 
where the new term / depends also on the fit coefficients 
c oi Fq. In what follows, we consider a single nucleus 
with the ground-state energy 



(34) 



(35) 



Thus, the new term is perturbatively small (assuming 
that the new fit coefficient cj is of "natural" order one). 
Furthermore, we make the assumption that / is of "nat- 
ural" size and perturbatively small in a neighborhood of 
{n*,z*) (i.e., its derivatives in this neighborhood are of 
the same order as /). We consider F — Fq + Cnew/new, 
and it is clear that the minimum of F will be found at 
some new occupations n* + 5n, z* -\- Sz with corrections 
Sn and Sz that are much smaller than n* and z* , respec- 
tively. We expand 



min [Fo{c;n,z) + Cff{c;n,z) 



Fo{c; n* + Sn, z* + Sz) + Cff{c; n* + Sn, z* + Sz) 
Fo{c; n*,z*) +Cff{c; n*,z*) 



dn 



K 

dnr 



Snp + ^ 



dFo 

dz„ 



Snp + CfY^ 



dZjy 



Sz^ 



Fo{c; n*,z*) +Cff{c; n*,z*), 



(36) 



and the approximation is due to our limitation to first- 
order corrections. Note that we have expanded to first- 
order in smallness, eliminating terms that go as deriva- 
tives of / since \fSn\ < |Fo| and \fSp\ < |i^o|- The 
derivatives of the functional Fq with respect to the occu- 
pation numbers vanish at the optimum occupation num- 
bers (n*, z*), where it is understood that the variation is 
only with respect to those occupation numbers that are 
not at any of the boundaries (i.e., those occupation num- 
bers for which n* ^ or dj and z* ^ or dj), and that 
the variation fulfills the equality constraints. Technically 
speaking, these are the reduced derivatives [IS]. Thus, 
in leading order of perturbation theory, the functional is 
simply a mass formula (as it is evaluated at the leading- 
order occupations), and the chi-square fit cannot yield 
an increased root-mean-square error for the ground-state 
energies. In the worst case, c/ = will result from the 
fit. Thus, only the addition of terms to the functional 
that cause perturbatively small changes to the occupa- 
tion numbers are expected to result in a decreased chi- 
square. 

For quickly evaluating a new candidate term / to be 
added to the functional, we begin by treating the full 
functional as a mass formula. That is, we freeze the 
proton and neutron occupations at the optimal values 
determined prior to the addition of the new term and 
approximate the ground-state energy in the presence of 
the new term as 

Ei{c;n*,z*) sa F{c;n* , z*) + Cf f^^^^in* , z*). (37) 



We evaluate the energy (37 1 for each nucleus i, perform a 
minimization of the resulting mass formula ( 37 1 , and 
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obtain a new set of test coefficients, ctcst- Using these 



test ' ^test 



and 



coefficients we recalculate the occupations Ti- 
the ground-state energies 

-E-ilctcst; 'T-tcsti ^tcst) ~ 

min (F{ctcsun,z)+Cff{ctcsun,z)) (38) 

that minimize the functional. The computation of the 
corrections {Sn, Sz) « («tost ~ -^tcst ~ ^*) to the occu- 
pation numbers shows whether the new term is pertur- 
bative in character and can thus be expected to lower the 
chi-square. 



E. Model for Nuclear Radii 

In DFT the radius r of a nucleus is computed from the 
density p(r) as 



i I dv' rV(r) 



(39) 



In constructing the functional we employed a harmonic 
oscillator basis for the shell-model and the scaling argu- 
ments. Consequently, we also employ it for the computa- 
tion of (r^). The expectation value of the radius squared 
in the harmonic oscillator shell with principal quantum 
number ly and angular momentum A is 



(40) 



Here, £ is the oscillator length and is set to 
V 492.5^1/'^ fm, and m is the mass of the proton. We 
use p = 2i> + X and thus find for the expectation value of 
the charge radius squared 



p 



(41) 



In computing charge radii, we must account for the finite 
size of the nucleons [46] , 



' ch,p/ 



N 



(42) 



Here, the proton and neutron charge radii are 



^\ ) = 

ch,pi 



0.877 fm. 



= —0.1161 fm^, respectively. 



We follow Ref. 47 and compute the charge radius 
within the following model: 



V = Vi + + v^- 



N - 



U4- 



(7V-Z)2 



(43) 
(44) 



Here r = VT^^c/iT' the coefficients v parameterize the 
model. We perform a fit of the radii to determine 
the coefficients v = (wi, 
whereby we minimize 



vq) in Eqs. (43) and (44 1, 



1 



^ 2—1 



Npts 

^ irfit4v; z) - r'T^'P f . 



(45) 



III. MINIMIZATION OF THE FUNCTIONAL 

The POUNDerS algorithm was developed for nonlin- 
ear least-squares problems where the derivatives of the 
residuals are unavailable. A summary of the algorithm 
in the context of DFT can be found in [3D] . 

For fixed n and z, the derivatives \/ cEf^{c;n, z) are 
known and continuous except when Cgs = —A^^^ for some 



nucleus; see Eq. (14). However, the form of the nonlinear 



lower-level minimization in Eq. ( 27 ) does not satisfy stan 



dard regularity conditions that would ensure existence 
and continuity of the derivatives ^-^-^ and '^^J^ (see, 

e.g., [IH]). Thus, unavailability of the residual derivatives 
in our case comes from the dependence of the optimal oc- 
cupation numbers n* , z* on the coefficients c. 

From a theoretical standpoint POUNDerS requires 
continuously differentiable residuals. For this functional, 
however, we found that a smooth model-based method 
accounting for the problem structure yielded better co- 
efficients in fewer simulations than do optimization al- 
gorithms that can explicitly treat the nonsmoothness. A 
similar result was found for many of the piecewise smooth 
problems in [49] . 



To minimize Eq. (29), the POUNDerS used in [30] was 



modified to account for known partial derivatives with 
respect to some of the coefficients. The optimal occu- 
pations {n*(c)}, {z*(c)} are independent of the three co- 
efficients {cc,cp^Ci) because the corresponding terms in 
the functional Eq. ( 13 ) do not involve n, z. Hence we can 



algebraically compute the (continuous) derivatives 



2 Rith 



dcj ' dcj dck ' 



j,fce {c,P,i}. 



(46) 



The unavailability of derivatives makes optimization 
significantly more challenging. Over the course of the 
optimization, POUNDerS effectively builds up coarse ap- 
proximations to first- and second-order derivatives of the 
residuals by interpolating the residuals at previously eval- 
uated coefficient values. Knowing the residual derivatives 
of 3 of the 17 coefficients effectively lowers the dimension 
of the difficult derivative-free optimization problem, re- 
sulting in fewer evaluations of x^- 

To guard against the effects of multiple local mini- 
mizers and discontinuities of the computed energies, we 
found that a sufficient strategy was periodic restarting 
of POUNDerS in neighborhoods of mild size. This al- 
lows the local algorithm to occasionally look beyond 
the smaller neighborhood it has focused on. While it 
is impossible to guarantee that has been globally 
minimized, the local solution reported in the next sec- 
tion is significantly better than the values found for 
other coefficients. Based on 50,000 c values uniformly 
drawn from the hypercube [—1,1]^'', Fig. [l shows that 
roughly 0.01% of c have < 10^ and roughly 50% have 
X^ > 10^. This shows that by taking into account the 
structure of and the availability of some derivatives. 
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fit coefRcients in units of MeV are 



10" 



lio- 



10" 



.« 10 
Q 



E 

LU 



10" 



10" 



/ 



10 



10 



10 



10 ■ 



FIG. 1: (Color online) log — log empirical cumulative distri- 
bution function showing the probability of randomly finding 
a value below the given value in the hypercube [—1, 1]^^. 



with POUNDerS we were able to find a proverbial "nee- 
dle in a haystack" with x ~ 1. 



IV. RESULTS 



We now analyze the results of the functional based on 
coefficients c determined by the ^ fit ( 29 1 of the energy 
calculations to experimental data. 



A. Energy 



We fit our functional to a set of 2,049 nuclei from the 
2003 atomic data evaluation [50] whose uncertainty in 
the binding energy is below 200 keV. The resulting fit 
produces a least-squares error of x = 1-31 MeV, and the 



Cc 
Cp 

Cs 
Cas 
^ss 
C2 
C3 
C4 

C5 
C6 
Cl 
Cs 
Cs 
C7 
C9 
ClO 
Cll 



-0.619948, 

11.170908, 

0.891816, 

7.434098, 

0.397623, 

23.174559, 

0.233345, 

0.493533, 

-10.678202, 

-0.353447, 

-7.672829, 

-0.429029, 

7.333364, 

0.112605, 

0.382828, 

-4.107110, 

1.383515. 



Our root-mean-square deviation is competitive with cur- 
rent Skryme force-based functionals [30]. While the 
current least-squares deviation is greater than the best 
achieved by a mean-field based functional , we provide 
a novel mass functional that uses the orbital occupations 
of nucleons as the relevant degrees of freedom. 

Figure [2] shows an iV vs. Z chart of the fit nuclei with 
color representing the difference i?"^ — E°^p between the 
energies computed from the functional and the experi- 
mental data for the 2,049 nuclei employed in the fit. The 
differences are a smooth function of the neutron and pro- 
ton numbers. Figures [3] and [4] display the energy differ- 
ences as functions of N and Z. There are still systematic 
deviations associated with shell oscillations. Recall that 
the shell closures are input to our functional through the 
choice of the single-particle degrees of freedom. Thus, 
the shell oscillations reflect smaller deficiencies associated 
with the description of nuclear deformation. 

To test the extrapolation properties of our functional, 
we fit the functional to a smaller set of 1,837 nuclei taken 
from the 1993 atomic evaluation data set [51]. For this 
data set, the root-mean-square error is x = 1-38 MeV. 
This least-squares deviation is close to that from a fit 
to the larger set of 2,049 nuclei, and the deviations are 
again smooth across the nuclear chart. Employing the 
functional from the fit to the 1993 data set, we compute 
the ground-state energies of all 2,049 nuclei and find a 
root-mean-square deviation of x = 1-34 MeV. Thus, the 
fimctional has good extrapolation properties. 

To further test the functional's predictive power, we 
use the coefficients c from the fit to the 1993 data set (x = 
1.38 MeV) and the coefficients from the fit to the 2,049 
nuclei (x = 1-31 MeV) and compute the binding energies 
of 2,149 nuclei in the complete 2003 nuclear data set. The 
functional fit to the 1993 data set yields x = 1-40 MeV, 





FIG. 2: (Color online) Chart of 2,049 well-measured nuclei 
from the 2003 atomic data evaluation in Ref. [5D], color show- 
ing the difference E^^ — E"^"^ between the calculated energy 
and experimental energy. The energy difference exhibits a 
smooth behavior across the whole chart. Some overbinding is 
seen in the area of very heavy nuclei and tin isotopes. Color 
ranges from -9.962 MeV (blue) to 4.340 MeV (red). 




120 



FIG. 4: Energy difference - as a function of Z for 

the same nuclei as in Fig. [2] The small oscillations around 
zero indicate a good description of nuclear shell structure. 




FIG. 3: Energy difference S*"" - E'^''^ as a function of N for 
the same nuclei as in Fig. |2] The small oscillations around 
zero indicate a good description of nuclear shell structure. 




FIG. 5: (Color online) Energy difference E*"^ - E'"'^ between 
theoretical and experimental results for the functional that 
is fit to the 1993 data set and applied to the subset B of 
well-measured nuclei of the 2003 data set. The deviations are 
smooth across the nuclear chart and consistent with Fig. [2] 
Color ranges from -8.530 MeV (blue) to 4.114 MeV (red). 



and the functional from a fit to the 2,049 nuclei yields x = 
1.38 MeV. These values are close to x = 1-37 MeV that 
results from a fit of the functional to the full 2003 data 
set. Thus, the extrapolation properties of the functional 
are quite good. Table |l] summarizes the details of how 
our functional extrapolates from data set to data set. 
Figure [5] shows the differences between theoretical and 
experimental ground-state energies across the chart of 
nuclei employed for the extrapolation. 




B. Radii 

We now turn to the results for nuclear charge radii. 
We fit our six-parameter model (44) for the charge radii 
to the experimental data of 772 nuclei from Ref. [55] and 
obtain a least-squares deviation of Xr = 0.047 fm. The 



FIG. 6: (Color onhne) Energy difference E*^ - £™p between 
theoretical and experimental results for the functional fit to 
the 1993 data set and applied to the 2003 data set. The 
deviations are smooth across the nuclear chart and consistent 
with Fig. [2] Color ranges from -9.263 MeV (blue) to 6.183 
MeV {redy 
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TABLE I: Root-mean-square deviations of binding energies 
resulting from a global fit of the functional and from extrap- 
olation to larger data sets. Data set A consists of all nuclei of 
the 1993 mass evaluation [51] , data set C consists of all nuclei 
of the 2003 mass evaluation |50| . and data set B is a subset 
of data set C consisting of well-measured nuclei whose uncer- 
tainty in the mass is below 200 keV. The number of nuclei in 
each data set is denoted by A^pts. 



Data 
Set 




X (MeV) 


Fit 


Extrapolation to 


Data Set B 


Data Set C 


A 


1837 


1.38 


1.34 


1.40 


B 


2049 


1.31 




1.38 


C 


2149 


1.37 







resulting fit coefficients (with units of fm) are 

vi = 11.004005, 

V2 = 2.000870, 

W3 = --3. 248108, 

t>4 = -0.279495, 

= 0.775617, 

vq = -0.557746. 

Note that the data set contains both spherical and de- 
formed nuclei. Figure [7] shows the difference between 
the calculated and experimental radii. For a comparison, 
note that Dufio and Zuker state a least-squares error of 
about Xr — 0.01 fm for the charge radii of spherical nu- 
clei [17]. Radii from Skyrme functionals exhibit a root- 
mean-square deviation of about 0.025 fm [15]. Figure [s] 
shows the difference between computed and experimental 
charge radii as a function of neutron number. The out- 
liers seen in Fig. [8] correspond to isotopic chain of Tb and 
the elements Rb, Sr and Zr with neutrons = 60 — 62. 



FIG. 7: (Color online) Difference between experimental and 
theoretical charge radii for the set of 772 nuclei from Ref. [52] . 
Color ranges from -0.172 fm (blue) to 0.155 fm (red). 

Let us also study the extrapolation properties of our 
mass model. We fit the model (44) to the charge radii on 
a subset of 494 nuclei (taken from Ref. [S2]) chosen near 
the valley of stability. This yields a least-squares error of 



cr 

x' 

I- 




FIG. 8: Chart of difference between calculated and experi- 
mental charge radii for 772 nuclei as a function of A''. Lines 
connect isotopes. 



Xr — 0.046 fm. When using this model to compute the 
charge radii on the full set of 772 charge radii, we find 
a slightly increased Xr — 0.048 fm, indicating that the 
extrapolation is successful. 

For each nucleus, the occupation numbers enter the 
computation of the charge radius. We thus need to un- 
derstand how our model for the charge radius depends on 
the functional employed for the computation of binding 
energies; in other words, its dependence on the coeffi- 
cients c. A sensitivity analysis of our fit to binding en- 
ergies provides us with a confidence interval for each of 
the coefficients c. We take a randomly chosen sample of 
five sets of coefficients {ci, . . . , C5} within the confidence 
interval and recompute the structure (i.e., the occupa- 
tion numbers) and the binding energies. The resulting 
least-squares deviations for binding energy range from 
X = 1.34 MeV to 1.78 MeV. Subsequently, we compute 
the charge radii for the nuclei of interest (without refit- 
ting the coefficients v of our model for the radii). If we 
refit the coefficients v of our mass model and adjust them 
to the change in the coefficients c of the energy functional, 
the least-squares error for the radii changes by at most 
4%. Thus, we find the surprising result that the model 
for radii is relatively independent of the functional's fit 
coefficients c. 
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V. SUMMARY 

We constructed an occupation number-based energy 
functional for the calculation of nuclear binding ener- 
gies across the nuclear chart. The relevant degrees of 
freedom for the functional are the proton and neutron 
orbital occupations in the shell-model. A global fit of 
a 17-paramctcr functional to nuclear masses yields a 
least-squares deviation of x = 1.31 MeV for the bind- 
ing energies, and a simple six-parameter model for the 
charge radii yields a root-mean-square deviation of Xr = 
0.047 fm. The fimctional has good extrapolation prop- 
erties, evident from the application of the functional fit 
to data from the 1993 atomic mass evaluation to the nu- 
clei of the 2003 atomic mass evaluation. The form of the 
functional is guided by scaling arguments and by results 
from analytically solvable models. Isospin and surface 
correction terms in the functional proved to be impor- 
tant. These terms were determined by a systematic in- 
vestigation of correlations among possible terms for the 
functional and by an analysis of their perturbative behav- 
ior. Additional terms, and therefore lower least-squares 
deviations, may be obtained through further use of this 



method, as well as greater investigation of the micro- 
scopic contributions of higher-order effects in surface and 
radial terms. 
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